
require(Rcpp)
output_version <- '2021-02-25'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
ntl.zones <- subset(ntl.zones,select=c(OBJECTID,GID_0))

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

##
## Roads
##
road_lk_chad_shp = paste0(wkdir,'/local/gis_data/018_transportation/roads_panel')
roads <- read_sf(dirname(road_lk_chad_shp),layer=basename(road_lk_chad_shp))
roads <- subset(roads,select=c(COUNTRY,R2014,R2008,R2003,R1998,R1996,R1993,R1991,R1990,R1988,R1986,
                               R1985,R1984,R1983,R1976,R1973,R1971,R1969,R1968,R1967,R1965))

st_crs(ntl.zones) <- st_crs(st_as_sf(roads))

# intersection
int = st_intersection(roads, ntl.zones)

# project coords for length calc
int.prj <- int %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs")
  
# find out about the length of each line segment
int.prj$len_m = st_length(int.prj)

## output
out <- as.data.frame(st_drop_geometry(int.prj))
head(out)
save.dta13(out, paste0(wkdir,"/proc_data/ROADLENGTH_",output_version,".dta"))
